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Introduction 

f There exists an altitude regime in the atmosphere that is within the continuum do- 
/ main, but wherein the conventional Navier-Stokes equations cease to be accurate. The 
j altitude limits for this so called continuum transition regime depend on vehicle size and 
!j speed. 1 Within this regime the thickness of the bow shock wave is no longer negligible 
!| when compared to the shock stand-off distance and the peak radiation intensity occurs 
j within the shock wave structure itself. For this reason it is no longer valid to treat the 

ji shock wave as a discontinuous jump and it becomes necessary to compute through the 

j‘i shock wave itself. To accurately calculate hypersonic flowfields, the governing equations 
; must be capable of yielding realistic profiles of flow variables throughout the structure of 
a hypersonic shock wave. 2 

The conventional form of the Navier-Stokes equations is restricted to flows with only 
small departures from translational equilibrium; it is for this reason they do not provide the 
capability to accurately predict hypersonic shock structure. Calculations in the continuum 
transition regime, therefore, require the use of governing equations other than Navier- 
Stokes. In this paper several alternatives to Navier-Stokes are discussed; first for the case 
of a monatomic gas and then for the case of a diatomic gas where rotational energy must 
be included. Results are presented for normal shock calculations with argon and nitrogen. 

Navier-Stokes Equations 


The study of flow with translational non-equilibrium is largely concerned with at- 
tempts to solve the Boltzmann equation 


|(n/) + c, ±(nf) + = { j t {n /)} 
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where q is the particle velocity, / is the velocity distribution function f(ci,Xi,t), and Xi 
is the external force per unit mass acting on a particle. 3 Unfortunately, this equation is 
difficult to solve for the general case and approximate methods are required. 

One such method, the Chapman-Enskog expansion, is an expansion of the distribution 
- function about translational equilibrium in terms^of the JCnudsen niunber (£ = A / L) as 
follows 

/ = /o(l + £01 + + ... )■ (2) 

The zeroth-order solution (£ — ♦ 0) corresponds to translational equilibrium and the govern- 
ing equations reduce to the Euler equations. In retaining first-order terms in £ we obtain 
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the Navier-Stokes equations, 
given, as 


The Navier-Stokes equations in one spatial dimension are 
^ + £ (f + f „ ) = 0 , (3) 


where 



F = 


pu 

pu 2 + p 
puho 
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Here, p is the mass density, p is the pressure, and u is the velocity in the x direction. 
The internal energy per unit mass is denoted as e, with the total energy per unit mass 
eo = e + u 2 / 2 and total enthalpy per unit mass ho = eo + p/p. The shear stress and heat 
flux terms for this one spatial dimension are given as 
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where T is the translational temperature. 

It is important to note that the derivation of the Navier-Stokes equations from Chapman- 
Enskog theory proves, rather than assumes, that Stokes hypothesis (pa = 0) holds for a 
monatomic gas. 3 Furthermore, Chapman-Enskog theory allows for a theoretical means of 
determining viscosity. For an inverse power molecule with an intermolecular repulsive force 
given by 

„ K 
F = 


a relationship exists between the coefficient of viscosity p and the translational temperature 
and is given as 4 


T j 1 “ 
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Within the same set of assumptions, the coefficient of thermal conductivity may be written 
as 


k = 
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Burnett Equations 

Retention of the second-order terms in £ yields what are known as the Burnett equa- 
tions. These equations should produce more accurate solutions for flows with large devia- 
tions from translational equilibrium and-will allow meaningful calculations imthe contin- 
uum transition regime. The Burnett equations for one spatial dimension are of the same 
form as the Navier-Stokes equations except that the closure relations are given as 

T=Ti+T2 (10) 
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and the u>i 's and 9{ 's are constants determined from Chapman- Enskog theory and depend 
on the constant 77 in the inverse power model. 4 

One drawback of the Burnett equations, other than their complexity, is that they 
tend to become unstable as mesh spacing is refined. This problem has been overcome by 
adding selected terms from the third-order Super-Burnett equations that have a stabilizing 
effect. 5 Hence, the Burnett equations are solved as 
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and 
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Here, t 3 and q 3 are the selected higher-order Super Burnett terms added for stability. 
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Rotational Energy 

The Navier-Stokes and Burnett equations are valid for monatomic gases only. For 
diatomic gases some means must be included to account for molecular rotation. Inclusion 
of molecular rotation from first principles requires a distribution function involving not 
only the coordinates of the molecule and its linear velocity, but also its rotational energy. 

The corresponding Boltzmann equation must also be modified so that the right-hand side 
includes a mechanism for energy exchange between-the translational and rotational energy- 
modes during collisions. 4 Wang Chang and Uhlenbeck 6 were the first to carry out the 
Chapman-Enskog expansion to first order for such a system. For cases where the relaxation 
time for rotation is sufficiently small such that rotational and translational energies axe at 
all times near their equilibrium values, the system yields the Navier-Stokes equations plus 
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an additional bulk viscosity term (jubV • V). This case is termed easy exchange and for 
one spatial dimension results in 

t = t 1 +t ib , (18) 


where 
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and q = qi, where 
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For cases where the relaxation time for rotation is significant and large deviations from 
rotational equilibrium occur (termed difficult exchange ) the system results in the Navier- 
Stokes equations with pe o = p(e lr + e r + u 2 / 2) plus an additional equation describing the 
relaxation of rotational energy 


de T e r — e r 
dt Z r T c 


( 21 ) 


Here, e r is the rotational energy at thermal equilibrium. The vibrational energy contribu- 
tion to the internal energy has been neglected for two reasons. Firstly, only cases where 
the initial conditions are such that vibrational energy is negligible are considered, and 
secondly, for most shocks, the vibrational energy mode requires many collisions to reach 
equilibrium (on the order of 10°) and is essentially frozen throughout the shock. Electronic 
energy contribution to the internal energy has also been neglected. The value Z r is known 
as the rotational collision number and the functional form is given as 11 
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and r c is the mean collision time given as 



(23) 


For one spatial dimension the governing equations for difficult exchange are given as 
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Here, r = T\, q = qi, where 

k = —pCvtf (26) 
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where 


k r = ncy r . (28) 

At present no one has carried out the Chapman-Enskog expansion to second order for 
a gas possessing rotational energy. 4 However, one obtains reasonably good results if the 
Burnett equations are used in conjunction with either of the two cases above for easy or 
difficult energy exchange. 


Bulk Viscosity Model 


As previously mentioned, when the energy exchange between rotation and translation 
occurs freely (r r is sufficiently small) the presence of rotational energy leads to the addition 
of a bulk viscosity term. This term appears in the stress components of the Navier-Stokes 
equations as follows 
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where p is the pressure corresponding to thermal equilibrium at T. It is important to note 
that bulk viscosity becomes important only when compressibility effects are present since 
/ic appears only with the divergence of the velocity vector. 3 From kinetic theory we have 

Oij = - < c x cj > , (30) 


and 


P = 3 P < Ci > ’ 


6tr — 2 < ^ ' 


(31) 


(32) 


Here, the brackets <> represent a moment of the distribution function, c { is the thermal 
velocity, and p is the kinetic pressure which is valid in thermal non-equilibrium. It follows 
that 

(33) 


and likewise 
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where e tr is the translational energy at'thermal equilibrium. It mow becomes apparent 
that the bulk viscosity is a relaxation phenomenon resulting from the lag in adjustment to 
equilibrium of the rotational mode. 7 For example, when a gas is expanded there will be a 
reduction in the internal energy (e). This reduction will initially come entirely from the 
translational modes (fast modes) and only after some relaxation time will the translational 


and rotational modes reach equilibrium at T. During this time the kinetic pressure is lower 
than that for thermal equilibrium and this is consistent with the equation above since 
(V ■ V) > 0. A similar process occurs for a compression. 

To develop an expression for the bulk viscosity, consider a small volume of gas (8Vol ) 
undergoing an expansion or compression as in Ref. 7. Considering only the work done by 
pressure the change in internal energy in a short time t* is given as 

-£(V-V>*. (35) 

P 

It is noted that in thermal equilibrium, for cases where vibrational and electronic energy 
contributions to the internal energy are negligible, e tr = 3/5e and e r = 2/5e. In an 
expansion where thermal equilibrium is maintained e tT decreases by the amount 

-ll(y.V)t* (36) 

5 P 

and e r decreases by the amount 

-7 -(V-V)t*. (37) 

5p 

In an expansion where there is a lag in the rotational mode the entire loss of energy will 
initially come from the translational mode such that 

e t r = e lr -l^(V-V)t\ (38) 

5 p 

In comparison with equation (34) in terms of bulk viscosity it follows that 

Ps = — P t • (39) 

It can be shown 7 that t* is approximately the relaxation time for rotation and taken as 

r r = Z r T c , (40) 

which is the relaxation time in the rotational model of the previous section. 

Numerical Issues 

_ An implicit Steger- Warming type flux-splitting algorithm as -in Ref. 8 was used to 
march in time to the steady-state solution. In Ref. 8 a general splitting was developed 
for a one-dimensional, vibrationally-relaxing, chemically reacting flow where it was shown 
that 

F = AQ = [T(A + + A - )T -1 ]Q = F + + F~ . (41) 
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Here, Q is the vector of conserved variables, and F is the flux vector given as 
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A is the flux Jacobian dF/dQ, A is the diagonal matrix whose elements are the eigenvalues 
of A, and it was shown that 


F = (— z — )p^A F a + ~p\d~Fb + —p^c F c, 


where 


F, = 


7 

/ Pi/P \ 
P2/P 

Pn/p 
u 

Pl e nll P 

PM^tiM / P 
\h 0 - a 2 /( 7 - 1 )/ 


27' 


F B,C = 


27' 

/ pl/p \ 
P2/P 

Pn/p 
udc a 

Pl^-nl/ P 

Pm^tim/p 
\ ho ±ua ) 


(43) 


(44a, b ) 


and Xa = u, \b = u + a, and Ac = u — a. For the present application the rotational 
non-equilibrium model is used in place of the vibrational model of Ref. 8. 

Both the Navier-Stokes and Burnett viscous fluxes are treated implicitly in a method 
similar to that of Ref. 9. Here, the coefficients in front of the derivative terms are frozen 
during the linearization. The source terms resulting from non-equilibrium rotation are also 
handled implicitly. An Euler implicit time integration method has been applied to the set 
of governing equations, resulting in the following semi-discrete form 
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and 


R(Q n ) = S X (F + + F“ + F v )” - W a n Ax. 


(48) 


A third-order, upwind biased extrapolation (in MUSCL form) is used for the inviscid fluxes 
in the right-hand side of equation (45), while a first order upwind extrapolation used for 
the inviscid Jacobians appearing in the left-hand side of equation (45). No limiting was 
necessary for any of the calculations presented. It is noted that Navier-Stokes calculations 
require a block tridiagonal inversion at each time step while Burnett calculations require 
a block pentadiagonal inversion at each time step. 

Each solution was obtained on a grid consisting of 150 equally spaced mesh points. 
The initial condition was composed of three segments, a constant region initialized at the 
inflow conditions, an initial guess at the shock profile that blended between the inflow 
conditions and the conditions behind a normal shock, and finally a constant region at the 
conditions behind a normal shock. Characteristic boundary conditions were used at the 
outflow with the pressure at the exit set at that behind a normal shock. 

In all the calculations presented the following expression is used as a reference length 
for shock structure 


16 /xi 

5pi(27r/?T 1 ) 0 - 5 ' 


(49) 


This is the mean free path that would exist upstream of the shock if the gas were composed 
of hard elastic spheres. Also, for all calculations the coefficients in the Burnett equations 
are those corresponding to a hard sphere gas (rj — oo,uj = .5) while u> = .72 is used in 
the viscosity law for better agreement with experiment. In nitrogen calculations where the 
rotational non-equilibrium model is incorporated the constants = 18. and T* = 91.5AT 
are used. 4 

Density and temperature profiles are presented in a normalized fashion where 


Pn — 


P~ P2 
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Results 

Results for argon shocks at Mach 20 and 35 are shown in Figures 1 and 2. Here, 
normalized temperature and density profiles are presented for Navier-Stokes and Burnett 
calculations as well as a Direct Simulation Monte Carlo solution (DSMC) as shown in Refs. 
1 and 9. It is assumed that for these conditions the DSMC solution is essentially the valid 
solution. The density thicknesses predicted by the Burnett- solution eompare-well with- the 
DSMC results in both the Mach 20 an 35 cases. The Navier-Stokes density thicknesses are 
slightly less than the DSMC results. In the case of temperature thicknesses, the Burnett 
solution represents a significant improvement over the Navier-Stokes solution. Overall, it 
is seen for both cases that inclusion of the Burnett terms yields a significant improvement 
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over Navier-Stokes for both temperature an density profiles. Our results compare well with 
the shock calculations from Refs. 1 and 9. 

Results for diatomic nitrogen at Mach 11 are shown in Figures 3 and 4. Figure 3 
shows calculations with the Navier-Stokes equations alone, the Navier-Stokes equations 
with rotational non-equilibrium, and the Navier-Stokes equations with bulk viscosity, as 
well as the DSMC result from Ref. 4. Again, the DSMC solution is assumed to be 
valid. It is seen that the Navier-Stokes calculation alone produces density thicknesses 
that are too small. Navier-Stokes calculations with bulk viscosity and rotational non- 
equilibrium produce density profiles that compare very well with the DSMC calculation. 
In the case of temperature thicknesses, the Navier-Stokes calculation with bulk viscosity 
represents a slight improvement over Navier-Stokes calculations alone. The calculation 
with rotational non-equilibrium is able to capture the peak temperature in the shock, 
however, it is noted that in all cases the Navier-Stokes terms do not predict large enough 
temperature thicknesses. Ref. 4 did not present Navier-Stokes results for this case. 

Figure 4 shows calculations with the Burnett equations alone, the Burnett equations 
with rotational non-equlibrium, and the Burnett equations with bulk viscosity, as well as 
the DSMC calculation from Ref. 4. It is seen that the Burnett calculation alone produces 
density thicknesses that are too small. Burnett calculations with bulk viscosity produce 
density profiles that compare very close to the DSMC calculation, and appear to be slightly 
better than the calculation with rotational non-equilibrium. 

In the case of temperature thicknesses, the Burnett calculation with bulk viscosity 
represents a slight improvement over Burnett alone, and the calculation with rotational 
non-equilibrium very closely predicts the value and location of the peak temperature in the 
shock. The Burnett calculation with rotational non-equilibrium most closely matches the 
DSMC results. Our results for this case agree with the results from Ref. 4. Figures 5 and 6 
show similar results in a Mach 6 shock for Nitrogen and these calculations agree well with 
the results from Ref. 10. In each case the Burnett calculations arc more accurate than 
the Navier-Stokes calculations. The results obtained using the Navier-Stokes equations or 
Burnett equations in conjunction with a bulk viscosity model represents a new contribution. 
It is expected that the bulk viscosity results would be in better agreement with DSMC 
solutions for cases with smaller amounts of rotational non-equilibrium. 

Conclusions 

In all cases the Burnett terms represent a considerable improvement over Navier- 
Stokes for both density and temperature profiles. For cases with rotational energy the 
Burnett equations with a bulk viscosity model represent an improvement over the Burnett 
equations alone, but do not contain a mechanism to predict the peak temperature within 
the shockwave. Navier-Stokes calculations with rotational non-equilibrium can predict 
the value and location of the peak temperature in the shock but suffer from inaccurate 
shock thickness. Burnett Calculations with rotational non-equilibrium produce the most 
accurate results of the continuum models tested. 
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